Selección de sitios para el establecimiento de una red de estaciones hidrométricas usando decisión multicriterio y análisis hidrológico con modelo digital de elevaciones. Suplemento metodológico
José Martínez
Michela Izzo
Introducción
A pesar de la importancia de los datos hidrométricos para diseñar e implementar políticas efectivas de gestión de recursos hídricos, actualmente, la red de estaciones hidrométricas en funcionamiento en República Dominicana es limitada y enfrenta diversos desafíos (Burn 1997; INDRHI 2019; A. K. Mishra y Coulibaly 2010; Ashok K. Mishra y Coulibaly 2009). La distribución de estas estaciones es poco uniforme y, las que se encuentran en funcionamiento, están amenazadas por problemas de mantenimiento y reposición, debido a los bajos presupuestos asignados para este propósito y a la falta de personal estable (INDRHI 2019).
En el año 2019, solo 30 de un total de 170 estaciones hidrométricas inventariadas se encontraban en funcionamiento (INDRHI 2019). Se desconoce el estado actual de la red, pero a la fecha de elaboración de este informe, el personal técnico del INDRHI indicó que algunas estaciones podrían haber salido del sistema por problemas de nombramiento de personal y mantenimiento. Por lo tanto, es esencial expandir esta red para mejorar la representatividad y la precisión de los datos hidrométricos, lo que constituye una prioridad crítica.
La selección de sitios adecuados para la instalación de estaciones hidrométricas es crucial para la obtención de datos precisos y representativos (Ashok K. Mishra y Coulibaly 2009). Se deben considerar varios criterios relevantes, como la heterogeneidad climática y topográfica de la región, ya que estas variables influyen en la distribución de las precipitaciones y, por extensión, en la escorrentía dentro de la cuenca (World Meteorological Organization (WMO) y The International Association of Hydrological Sciences 1976). Además, la accesibilidad del sitio es fundamental para garantizar que el mantenimiento y la operación de las estaciones se pueda llevar a cabo de manera eficiente (Rojas Briceño et al. 2021).
Los criterios arriba mencionados fueron usados por este mismo equipo de investigación para proponer una red de estaciones meteoclimáticas, por lo consideramos oportuno aprovecharlos igualmente para diseñar la red de estaciones hidrométricas. En particular, las estaciones hidrométricas deben colocarse sobre corrientes fluviales donde midan los caudales de agua, por lo que es necesario contar con información geográfica precisa y densa sobre la red hidrográfica dominicana. Varios criterios adicionales deben cumplir los sitios de medición (Rantz 1982a, 1982b):
- El tramo de corriente o curso fluvial (río, arroyo, cañada) elegido debe ser recto durante unos 100 metros corriente arriba y abajo del sitio de medición.
- En todos los niveles, el flujo total está en una única corriente, sin flujo subterráneo que evite el sitio de interés.
- El lecho del curso fluvial y sus márgenes, no sufren erosión ni sedimentación y están libres de plantas acuáticas (sección estable).
- Las orillas de la corriente fluvial son permanentes, sin matorrales y lo suficientemente altas como para contener inundaciones.
- Preferiblemente, el curso tiene controles naturales invariables (sección estable): son preferidos los lugares con afloramientos de sustrato y con rápidos estables para condiciones de flujo bajo, así como los estrechamientos del canal o una cascada o caída que no esté sumergida en momentos de flujo alto. Cuando estas condiciones no se pueden conseguir en el sitio elegido, entonces es necesaria la realización de una obra, como por ejemplo un vertedero. Sin embargo, no en todos los ríos es posible realizar esta obra, sobre todo si se trata de ríos muy anchos, por lo que se recurre a métodos alternos de medición de velocidades y estimación de área de la sección (área elemental).
- En niveles extremadamente bajos, existe escorrentía o pozas aunque sólo sea aguas arriba, lo cual permite la continuidad de medición y evita altas velocidades en momentos de restablecimiento del caudal.
- El sitio candidato está lo suficientemente alejado de confluencias o de los efectos de las mareas para evitar impactos indeseados en la medición.
- Hay un tramo disponible para medir el caudal en todos los niveles cerca del sitio de medición.
- El sitio es accesible para instalar, operar y mantener el lugar de medición, siendo su elección un compromiso entre estos criterios.
Para obtener algunos de estos criterios, es necesario compilar datos hidrológicos, geológicos y topográficos a nivel nacional. En principio, y reduciendo el análisis de selección de sitios a fuentes geoespaciales disponibles o derivables, se deben identificar tramos de cursos fluviales que cumplan con las siguientes características (Rantz 1982a, 1982b):
- Alineación recta.
- Rocas consolidadas expuestas (en contraposición a los sedimentos aluviales).
- Las márgenes sujetas a desbordamiento.
- Márgenes pronunciadas para garantizar el confinamiento del flujo.
- Canales divididos.
- Posible impacto de fenómenos/estructuras aguas arriba que induzcan variabilidad. Estos efectos pueden ser causados por afluentes, confluencias o embalses.
- Viabilidad de instalar medidores de flujo.
Los criterios anteriores pueden derivarse a partir de una precisa y densa red hidrográfica del país, además de disponiendo de datos complementarios que puedan manejarse de forma integrada por medio de un sistema de información geográfica. Sin embargo, las fuentes de información geográfica sobre la hidrografía dominicana disponibles al público, no cuentan con la resolución espacial requerida como para generar todos los criterios relacionados arriba.
Existen al menos dos fuentes de información geográfica sobre la red hidrográfica dominicana. La primera es la red digitalizada a partir del mapa topográfico nacional escala 1:50,000 (“MTN-50k”) (Instituto Cartográfico Militar (ICM) 1989). La red del MTN-50K, aunque es bastante exhaustiva, es realmente una red parcial e inconsistente (e.g. uso de criterios subjetivos para definir cursos fluviales), y no aporta características morfométricas de las corrientes (e.g. jerarquía), además de que no es lo suficientemente precisa como para utilizarla, por sí sola, en análisis hidrológicos. No obstante, dicha red puede servir de insumo para crear el DEM con fines hidrológicos, como por ejemplo, la aplicación o tallado de la red.
La segunda fuente sobre la hidrografía dominicana es un conjunto de documentos técnicos y multitemáticos, tanto de ámbito nacional (INDRHI 2012, 1996; INDRHI y AQUATER 2000; INDRHI y EPTISA 2000; OEA y INDRHI 1994; Secretaría de Estado de Medio Ambiente y Recursos Naturales 2004; Rodríguez y Febrillet 2006), como de cuencas seleccionadas (CIDIAT y INDRHI 1992; Halcrow-COR Ing. S.A. 2002; SERCITEC y INDRHI 2002), realizados mayormente por personal del Instituto Nacional de Recursos Hidráulicos (INDRHI). Sin embargo, cada estudio fue elaborado con métodos distintos, por lo que no es viable la consolidación de una hidrografía consistente de ámbito nacional a partir de estos. Cabe destacar que se están dando pasos en esa dirección actualmente, tras haberse iniciado la actualización del Plan Hidrológico Nacional.
El presente estudio tiene, como primer objetivo generar una red densa de ríos, arroyos y cañadas (en la terminología geomorfológica, talweg, o línea que une puntos de elevaciones mínimas en vaguadas), elaborada a partir de las fuentes geoespaciales abiertas más precisas de las que se tenga conocimiento, con miras a seleccionar sitios idóneos para la instalación de estaciones hidrométricas. Para este cometido, optamos por generar la red a partir del modelo digital de elevaciones de mayor resolución espacial disponible al público en la actualidad, usando métodos semiautomáticos implementados en software de código abierto, y apoyándonos en fuentes pre-existentes, como el mapa topográfico nacional y el mapa geológico, entre otros.
Como segundo objetivo nos planteamos seleccionar sitios idóneos para la instalación de estaciones hidrométricas, considerando únicamente la escorrentía superficial en talwegs, excluyendo mediciones de pozos, lagos, canales de riego, embalses y otros cuerpos de agua distintos a ríos, arroyos y cañadas. Para esta labor, y como parte del diseño de la investigación, aprovechamos las unidades hexagonales ya elegidas previamente en el estudio complementario a este, donde se propuso aumentar la densidad de la red estaciones meteoclimáticas (ver informe titulado “Selección de sitios para el establecimiento de una red de estaciones meteoclimáticas en República Dominicana usando decisión multicriterio y análisis de vecindad”).
Notas sobre la unidad elemental talweg y la precisión de la red. Un talweg es una línea imaginaria sobre el terreno que une los puntos de más bajos de un valle (Foucault y Raoult 1985). Por lo tanto, no es el equivalente a un curso fluvial permanente, pues la circulación de agua dependerá de otros factores además de los meramente topgráficos (clima, umbral de acumulación elegido, sustrato, entre otros). En nuestro mapa, un talweg refleja el espacio por donde, potencialmente, el sobrante de agua circularía en forma de escorrentía superficial en condiciones de suelo saturado. Sin embargo, el trazado final que describa la escorrentía, estará condicionado por múltiples factores (temporales, de cobertura, de suelo) siendo en gran medida dependiente de la precisión del DEM y del relieve en cuestión. Particularmente, sobre la precisión de la red, consideramos que es alta en sistemas montañososos no kársticos, y media en sistemas kársticos. En particular, refiriéndonos a los sistemas kársticos, la incertidumbre fue ligeramente mayor en nuestro mapa porque no se dispone de un inventario exhaustivo y preciso de depresiones (e.g. dolinas), que son, en última instancia, las que condicionan la topología y jerarquización de la red en estos relieves. Para llenar este vacío, a partir del modelo de elevaciones elegido, creamos nuestro propio inventario de depresiones, pero lo usamos con cautela para definir la localización de sumideros reales por donde el flujo se infiltra al karst. Finalmente, es importante señalar que, en áreas urbanas y llanas, así como en campos con canales de riego, nuestra red puede tener una baja precisión. No obstante, esto no afecta el alcance de nuestro estudio, ya que nuestros objetivos se centran en otras áreas.
Notas sobre cuencas transfronterizas. Las cuencas transfronterizas fueron caracterizadas sólo con la información de la porción de territorio dominicano, por lo cual el orden de red obtenido para las redes correspondientes, no necesariamente refleja el orden de real. Estas cuencas incluyen las de los ríos Artibonito, Pedernales y Masacre, así como varias cuencas pequeñas distribuidas a ambos lados de la línea fronteriza. Sin embargo, una rápida inspección visual sugiere que la referida omisión, no afecta el objetivo de nuestro estudio.
Los resultados de este estudio ofrecen la oportunidad de mejorar la red de estaciones hidrométricas, y tienen el potencial de informar medidas de conservación, planificación y gestión de recursos hídricos a nivel nacional, en un momento en el que la humanidad y, en particular, la sociedad dominicana, se preparan para afrontar los efectos del cambio climático y la consecuente escasez de agua pronosticada.
Materiales, métodos y resultados
Datos, paquetes y funciones
Los siguientes bloques de código cargan los paquetes, funciones y datos necesarios para preparar el índice de escenas del modelo digital de elevaciones ALOS PALSAR, y para representarlas en el mapa.
source('R/funciones.R')
library(raster)
library(sf)
library(kableExtra)
library(tidyverse)
library(gdalUtils)
estilo_kable <- function(df, titulo = '', cubre_anchura = T) {
df %>% kable(format = 'html', escape = F, booktabs = T, digits = 2, caption = titulo) %>%
kable_styling(bootstrap_options = c("hover", "condensed"), full_width = cubre_anchura)
}
dem_proc_dir <- 'alos-palsar-dem-rd/dem/'DEM: obtención y preparación
Selección de los datos de elevación
Elegimos el Modelo Digital de Elevación del Satélite de Observación Avanzada de la Tierra (ALOS-DEM, por sus siglas en inglés), pues es actualmente la fuente abierta de mayor resolución espacial (12.5 m) actualmente disponible al público (ASF DAAC 2015). A pesar de su alta resolución espacial, investigadores encontraron que no necesariamente es el de mayor precisión en cuanto a los estadísticos básicos de elevación (e.g. RMSE comparado con estaciones GNSS) (Aziz y Rashwan 2022). No obstante, la evaluación realizada por este equipo de investigadores no realizó preprocesamiento de los datos descargado. Sin embargo, para la extracción de redes de drenaje, a mejor resolución, mejores resultados en cuanto a densidad de la red y caracterización morfométrica. Realizamos la descarga desde el Centro de Archivo Activo Distribuido (DAAC) del Alaska Satellite Facility (ASF). Usando este modelo de elevaciones como fuente, extrajimos de manera automática las redes hidrográficas de todo el país.
Descarga de escenas
Identificamos las escenas necesarias para cubrir el país usando una búsqueda geográfica mediante polígono delimitador. Dado que la misión del ALOS-DEM ofrece escenas de distinta fecha para una misma área, descargamos escenas redundantes que posteriormente excluimos del análisis. La descarga la realizamos por lotes, usando un script de Python provisto por el propio ASF.
python download-all-2023-04-20_00-30-00.pyCreación de índice de DEM
Utilizando el índice de huellas de escenas, escribimos un pequeño programa para seleccionar la escena más reciente en las áreas donde contábamos con escenas redundantes. Con esto construimos un índice de DEM seleccionados.
ind_orig <- invisible(
st_read('alos-palsar-dem-rd/asf-datapool-results-2023-04-19_08-31-26.geojson', quiet = T)) %>%
rownames_to_column('fila') %>% mutate(fila = as.integer(fila))
distancias <- ind_orig %>% st_centroid() %>% st_distance() %>% units::drop_units()
distancias[upper.tri(distancias, diag = T)] <- NA
indices <- which(distancias < 1000, arr.ind = TRUE)
duplicados <- as.data.frame(indices) %>%
mutate(dup_id = 1:nrow(indices)) %>%
pivot_longer(-dup_id, names_to = 'tipo', values_to = 'fila') %>%
select(-tipo)
seleccionados <- duplicados %>%
inner_join(ind_orig %>% select(fila, startTime) %>% st_drop_geometry) %>%
group_by(dup_id) %>% filter(startTime == max(startTime)) %>% pull(fila)
ind_orig_sel <- ind_orig %>% filter(!fila %in% duplicados$fila | fila %in% seleccionados) %>%
filter(centerLon < -72.1821)En total, para cubrir el territorio de República Dominicano, necesitamos 28 de escenas únicas ALOS-PALSAR.
ind_orig_sel %>% select(sceneName, startTime) %>% st_drop_geometry() %>%
kable(format = 'html', escape = F, booktabs = T,
caption = 'Escenas ALOS-PALSAR usadas para generar un DEM de 12.5 m de resolución espacial de República Dominicana') %>%
kable_styling(bootstrap_options = c("hover", "condensed"), full_width = T)| sceneName | startTime |
|---|---|
| ALPSRP253240380 | 2010-10-25 23:18:16 |
| ALPSRP253240370 | 2010-10-25 23:18:08 |
| ALPSRP253240360 | 2010-10-25 23:17:59 |
| ALPSRP253240350 | 2010-10-25 23:17:51 |
| ALPSRP252510370 | 2010-10-20 23:11:46 |
| ALPSRP252510360 | 2010-10-20 23:11:38 |
| ALPSRP252510350 | 2010-10-20 23:11:30 |
| ALPSRP251490380 | 2010-10-13 23:22:45 |
| ALPSRP251490370 | 2010-10-13 23:22:36 |
| ALPSRP251490360 | 2010-10-13 23:22:28 |
| ALPSRP251490350 | 2010-10-13 23:22:20 |
| ALPSRP251490340 | 2010-10-13 23:22:12 |
| ALPSRP250760380 | 2010-10-08 23:16:23 |
| ALPSRP250760370 | 2010-10-08 23:16:15 |
| ALPSRP250760360 | 2010-10-08 23:16:06 |
| ALPSRP250760350 | 2010-10-08 23:15:58 |
| ALPSRP250030360 | 2010-10-03 23:09:44 |
| ALPSRP250030350 | 2010-10-03 23:09:36 |
| ALPSRP248280370 | 2010-09-21 23:14:21 |
| ALPSRP248280360 | 2010-09-21 23:14:13 |
| ALPSRP248280350 | 2010-09-21 23:14:05 |
| ALPSRP247260360 | 2010-09-14 23:25:03 |
| ALPSRP247260350 | 2010-09-14 23:24:55 |
| ALPSRP247260340 | 2010-09-14 23:24:47 |
| ALPSRP242300380 | 2010-08-11 23:21:28 |
| ALPSRP242300370 | 2010-08-11 23:21:19 |
| ALPSRP242300360 | 2010-08-11 23:21:11 |
| ALPSRP242300350 | 2010-08-11 23:21:03 |
Representamos el índice de escena con ggplot2, superponiendo las huellas (polígono de área con datos) sobre el límite costero e internacional de República Dominicana.
FIGURA 1: Mapa índice de escenas usadas en la formación del DEM 12.5 ALOS-PALSAR de RD
Extracción de las escenas seleccionadas
Usando el índice de las escenas seleccionadas, extrajimos los correspondientes DEM en formato GTiff desde los archivos comprimidos (.zip).
zip_path <- 'alos-palsar-dem-rd/'
sapply(ind_orig_sel$fileName,
function(x)
unzip(
zipfile = paste0(zip_path, x),
exdir = paste0(zip_path, 'dem'), junkpaths = T,
files = paste0(gsub('.zip', '', x), '/', gsub('zip', 'dem.tif', x)))
)Verificación de CRS, transformar a zona 19N
Todos los DEM fueron servidos bajo el sistema de coordenadas universal transversal de Mercator (UTM), pero algunos fueron fijados en el huso 18N, por lo que los transformamos al huso 19N para generar un producto continuo usando la herramienta gdalwarp de la biblioteca GDAL (GDAL/OGR contributors 2022).
dems_orig_path <- list.files(path = 'alos-palsar-dem-rd/dem', pattern = '*dem.tif', full.names = T)
crs_18n <- names(which(sapply(dems_orig_path, function(x){
crs_x <- gdal_crs(x)
is_z18 <- grepl('zone 18N', crs_x[['wkt']])
})))
sapply(crs_18n, function(x) file.rename(from = x, to = gsub('.tif', '_z18n.tif', x)))
crs_18n_ren <- list.files(path = 'alos-palsar-dem-rd/dem', pattern = 'z18n.tif', full.names = T)
sapply(crs_18n_ren, function(x){
gdalwarp(
srcfile = x,
dstfile = gsub('_z18n.tif', '.tif', x),
t_srs = 'EPSG:32619', overwrite = T)})Construcción de ráster virtual
A efectos de eficientizar la manipulación del DEM, creamos un ráster virtual (VRT) usando la herramienta gdalbuildvrt de la biblioteca GDAL.
gdalbuildvrt(gdalfile = dems_orig_path,
output.vrt = paste0(paste0(zip_path, 'dem'), '/dem_seamless.vrt'),
resolution = 'highest', r = 'average')Creación de base de datos de GRASS GIS
Creamos una base de datos y localización de GRASS GIS usando como fuente de extensión y resolución el ráster virtual creado en el paso anterior (GRASS Development Team 2022). Decidimos usar GRASS GIS a partir de este punto para algunas tareas en las que este paquete es bastante eficiente (e.g. rellenado de nulos). Sin embargo, en pasos posteriores, alternamos el flujo de procesamiento con otras herramientas, como WhiteboxTools (John B. Lindsay 2018). En todo caso, nuestro criterio fue siempre aprovechar al máximo los recursos de hardware y software disponibles para obtener los productos requeridos en el menor tiempo posible.
# Usando Bash, desde la ruta ./alos-palsar-dem-rd/dem
grass --text -c dem_seamless.vrt ./grassdata
# Para abrir luego de cerrada: grass grassdata/PERMANENT/Máscara: crear, importar, aplicar
Creamos una máscara de país en QGIS (QGIS Development Team 2021), superponiendo el límite oficial obtenido desde la página de la Oficina Nacional de Estadística (ONE), y combinándolo con otras fuentes disponibles en línea, como GADM, Humanitarian Data Exchange (OCHA) y OpenStreetMap (Oficina Nacional de Estadística (ONE) 2018; GADM 2022; OCHA 2022; OpenStreetMap contributors 2017). De la máscara, eliminamos las superficies de máximas de lagos y lagunas no artificiales, pues nos interesa procesar las cuencas endorreicas que drenan hacia ellos. No obstante, los embalses no los incluimos en dicha superficie, dado que necesitamos construir la jerarquía de red ignorando su presencia, es decir, asumiendo como continuos todos los cursos fluviales. Sobre esta máscara, creamos un área de influencia, para recortar el DEM con un cierto “acolchado” que nos permitiera análizar sin dificultades las áreas costeras y de frontera.
La creación de esta máscara fue el único paso que realizamos de forma semimanual, pues el resto del flujo de procesamiento lo realizamos con algoritmo automáticos.
Importamos la máscara generada a la base de datos de GRASS y la aplicamos. GRASS opera de forma eficiente, aplicando los algoritmos sólo dentro del área definida como máscara. Las áreas fuera de ésta son excluidas, eficientizando los recursos y evitando malgastar tiempo de CPU en áreas que ajenas al proyecto.
# Importar máscara
v.import input=mascara-1km.gpkg output=mascara_1km
# Fijar máscara
r.mask -r
r.mask vector=mascara_1km
# Ver ambiente
g.gisenv
## GISDBASE=/media/jose/datos/alos-palsar-dem-rd/dem
## LOCATION_NAME=grassdata
## MAPSET=PERMANENT
## GUI=text
## PID=1632142Creación de mapa ráster a partir de VRT
Importamos el ráster virtual a la base de datos de GRASS GIS. Con este paso generamos un mapa ráster dentro de la base de datos GRASS GIS, el cual es una realización con celdas manipulables y a la que le podemos aplicar algoritmos ráster de nuestra preferencia.
# Importar DEM a región de GRASS
time r.import --overwrite input=dem_seamless.vrt output=dem
## real
# Ver en lista (q para salir)
g.list type=raster
FIGURA 2: DEM ALOS PALSAR sin procesar, representado como relieve sombreado. Nótesense los píxeles sin datos, destacados en color rojo (Los Patos-Ojeda-Paraíso, provincia Barahona, sudoeste de República Dominicana)
Rellenado de píxeles nulos (sin datos)
Para esta tarea, utilizamos el eficiente complemento de GRASS r.fill.nulls. Lo configuramos para rellenar píxeles nulos usando interpolación spline bilineal con regularización Tykhonov (spline es un método de descomposición de curvas en porciones descritas por polinomios).
# Rellenar vacíos
time r.fillnulls --overwrite --verbose \
input=dem method="bilinear" \
tension=40 smooth=0.1 edge=3 npmin=600 segmax=300 lambda=0.01 \
output=dem_relleno
# Enviar mensaje al finalizar (ejecutar conjuntamente con anterior)
echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 10m11.925s
FIGURA 3: DEM ALOS PALSAR sin procesar, representado como relieve sombreado. Los píxeles sin datos fueron eliminados (Los Patos-Ojeda-Paraíso, provincia Barahona, sudoeste de República Dominicana)
Suavizado preservando morfologías
Para el suavizado, usamos la herramienta FeaturePreservingSmoothing de WhiteboxTools, la cual reduce reducir la rugosidad generada por ruido en el modelo digital de elevaciones (John B. Lindsay, Francioni, y Cockburn 2019; John B. Lindsay 2018). Para aplicar esta herramienta, primero exportamos el DEM desde la base de datos de GRASS GIS a archivo GeoTIFF, y posteriormente aplicamos el suavizado.
# Exportar a GTiff con compresión LZW
time r.out.gdal --overwrite --verbose createopt="COMPRESS=LZW,BIGTIFF=YES" \
input=dem_relleno \
format=GTiff type=Float64 output=dem_relleno.tif
# Enviar mensaje al finalizar (ejecutar conjuntamente con anterior)
echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 0m58.924s
# Comenzó a 23.20 de 22 de abril
time ~/WhiteboxTools_linux_amd64/WBT/whitebox_tools \
--wd='/media/jose/datos/alos-palsar-dem-rd/dem/' \
--filter=25 --norm_diff=45 --num_iter=5 \
--run=FeaturePreservingSmoothing --input='dem_relleno.tif' \
--output='dem_relleno_suavizado.tif' -v
# Enviar mensaje al finalizar (ejecutar conjuntamente con anterior)
echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 9min46.103s
FIGURA 4: DEM ALOS PALSAR suavizado, representado como relieve sombreado. Nótese la conservación de las morfologías principales y la eliminación del ruido sobre éstas (Los Patos-Ojeda-Paraíso, provincia Barahona, sudoeste de República Dominicana)
Importar DEM suavizado a la base de datos de GRASS GIS
Importamos el DEM suavizado para aplicarle nuevos procesamientos hidrológicos.
time r.import input=dem_relleno_suavizado.tif output=dem_suavizado
echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 0m21.593sObtener alturas pseudo-ortométricas
Usamos el ráster de altura de geoide de La Española a 1 minuto de resolución (EGM2008) para obtener alturas pseudo-ortométricas, por medio de una suma algebraica simple. Sin embargo, previamente fue necesario aumentar la resolución del ráster de altural del geoide antes de realizar la suma. Para esto, usamos r.resamp.rst (evaluamos una alternativa 2 con el addon r.resamp.interp, pero resultó ): https://grass.osgeo.org/grass82/manuals/r.resamp.interp.html
# Importar DEM a región de GRASS
r.import --overwrite input=egm2008-1_espanola.tif output=egm2008_1min
# Ver en lista (q para salir)
g.list type=raster
# Ver atributos de la región
g.region -p
# Alternativa 1. Usando r.resamp.rst. Más eficiente y precisa
# Fijar la región al geoide importado
g.region raster=egm2008_1min -ap
# Realizar la interpolación
r.resamp.rst --overwrite input=egm2008_1min ew_res=50 ns_res=50 elevation=egm2008_hires
echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real
# Fijar región a nuevo geoide
g.region raster=egm2008_hires -ap
# Alternativa 2. Usando r.resamp.interp. También eficiente, pero eliminar áreas de borde
# g.region res=50 -ap
# r.resamp.interp --overwrite input=egm2008_1min \
# output=egm2008_hires method=bilinear
# Exportar para explorar visualmente
# r.out.gdal --overwrite --verbose createopt="COMPRESS=LZW" \
# input=egm2008_hires \
# format=GTiff type=Float64 output=egm2008_hires.tif
# Volver a resolución de DEM rellenado y suavizado
g.region raster=dem_suavizado -ap
# Aplicar álgebra de mapas
r.mapcalc --overwrite "dem_pseudo_ortometrico = dem_suavizado - egm2008_hires"
FIGURA 5: Alturas respecto de geoide EGM08 (~ortométrica) y sobre elipsoide WGS84, de un transecto descendente desde Bahoruco Oriental al Mar Caribe (Los Patos-Ojeda-Paraíso, provincia Barahona, sudoeste de República Dominicana)
Tallar (aplicar) la red sobre el DEM
El tallado, grabado o aplicación de red (stream burning) consiste en reforzar, sobre el DEM, la red hidrográfica conocida para garantizar que los algoritmos automáticos de análisis hidrológico conduzcan el flujo a través de lechos existentes. Este procedimiento es particularmente útil (e imprescindible) en áreas llanas, pues produce redes hidrográficas más reales y, en suma, mejora la topología de la red. No obstante, el procesamiento modifica sensiblemente el DEM, especialmente en los lugares donde ocurre el grabado.
Probamos tres alternativas de tallado de DEM usando las herramientas FillBurn de WhiteboxTools (John B. Lindsay 2018), r.carve de GRASS GIS y álgebra de mapas, también de GRASS GIS (GRASS Development Team 2022). Asimismo, con cada algoritmo probamos intentando grabar en el DEM dos redes de drenaje distintas, una densa y otra compuesta sólo los cursos largos.
Para generar la red densa, nos apoyamos en el mapa topográfico nacional a escala 1:50,000 (MTN-50K) (Instituto Cartográfico Militar (ICM) 1989), el cual se encuentra digitalizado pero procede fuente anónima (esperamos que este estudio nos ayude a establecer la autoría para incluirla oportunamente). En particular, la red densa consistió en una recopilación exhaustiva de los talweg dibujados en el mapa topográfico (lechos en valles drenados, líneas de mínima elevación, “vaguadas”), excluyendo canales de riego y otras obras de ingeniería de reconducción de flujo. Para garantizar la continuidad de la red, recuperamos el trazado de la red previo al llenado de embalses, usando mapas topográficos de la época pre-embalse. Importamos esta red a la base de datos de GRASS GIS (una versión topológicamente corregida con el complemento v.clean de GRASS GIS) y obtuvimos unos estadísticos básicos de la misma para fines de comparación con la red generada en este estudio.
# Importar red hidrográfica digitalizada desde el MTN-50K
v.import --overwrite input=red_mtn50k_cleaned.gpkg output=red_mtn50k_cleaned
# Ver mapa importado en lista (q para salir)
g.list type=vector
# Calcular y pasar a archivo, la longitud de cursos y número de segmentos (ejecutar en casos de actualización)
v.to.db -c -p option=length map=red_mtn50k_cleaned | grep total > stats_length_red_mtn50k_cleaned.txt
v.to.db -c -p option=count map=red_mtn50k_cleaned | grep total > stats_count_red_mtn50k_cleaned.txtA partir del archivo de texto generado, obtuvimos los estadísticos básicos: se trata de una red compuesta por 29750 segmentos que suman un total de 45698.8 kilómetros de longitud.
FIGURA 6: Mapa de la red densa de ríos, arroyos, cañadas y canales de riego, alegadamente extraída desde el MTN-50K. Lo que sabemos que es fue digitalizada a partir del MTN-50K. Sin embargo, aunque está bastante bien alineada en buena parte del país, en algunos puntos del país (e.g. costa norte, sierra de Bahoruco), la red se desplaza ligeramente. La utilizamos para tallar el DEM, pero los análisis hidrológicos produjeron talwegs duplicados y, en los ríos con meandros divagantes, el resultado no se adaptaba bien al DEM (la fecha del DEM es 2010, mientras que la red es del siglo pasado). Autoría desconocida (esperamos que nos la indiquen próximamente)
Por otra parte, tal como indicamos arriba, alternativamente probamos con una red de cursos largos, la cual generamos a partir del MTN-50K (Instituto Cartográfico Militar (ICM) 1989), imágenes satelitales, y (especialmente) mapas de distintos proveedores (Google; Airbus, CNES; Airbus, Landsat; Copernicus; Maxar Technologies; U.S. Geological Survey 2023; OpenStreetMap contributors 2017). Esta red consistió en una selección de los ríos más largos y represados de República Dominicana, e incluyó tramos de ríos que atraviesan áreas problemáticas, como grandes llanuras o relieves kársticos. Al igual que en la red densa, para garantizar la continuidad topológica de la red, ignoramos la presencia de los embalses y grabamos trazados históricos sobre el DEM. Importamos esta red a la base de datos de GRASS GIS y obtuvimos estadísticos básicos.
# Importar red a GRASS
v.import --overwrite input=red_mtn50k_cleaned_largos.gpkg output=red_mtn50k_cleaned_largos
# Ver mapa importado en lista (q para salir)
g.list type=vector
# Calcular y pasar a archivo, la longitud de cursos y número de segmentos (ejecutar en casos de actualización)
v.to.db -c -p option=length map=red_mtn50k_cleaned_largos | grep total > stats_length_red_mtn50k_cleaned_largos.txt
v.to.db -c -p option=count map=red_mtn50k_cleaned_largos | grep total > stats_count_red_mtn50k_cleaned_largos.txtA partir del archivo de texto generado, obtuvimos los estadísticos básicos: se trata de una red compuesta por 1047 segmentos que suman un total de 4895.876 kilómetros de longitud.
FIGURA 7: Mapa de la red de cursos largos creada para el estudio a partir de varias fuentes (más detalles, en el texto).
El siguiente consistió en probar las herramientas disponibles, r.carve de GRASS GIS, álgebra de mapas (r.mapcalc``de GRASS GIS) yFillBurn` de WBT, con cada una de las redes.
Alternativa 1. Tallado con GRASS GIS, complemento r.carve (descartada)
La herramienta r.carve de GRASS GIS fue diseñada para grabar el DEM sin modificarlo sensiblemente, permitiendo configurar la profundidad y la anchura del grabado. Por defecto, la anchura de lecho es equivalente a la resolución del DEM. La profundidad puede definirse por el usuario, para lo cual nosotros establecimos 100 metros. Al intentar tallar la red densa sobre el DEM con este método, dedicamos un tiempo de cómputo de más de siete horas (sin paralelización), por lo que nos vimos obligados a interrumpirlo sin concluir. No obstante, la red de cursos largos sí pudimos tallarla sobre el DEM con esta herramienta, generando un resultado óptimo, aunque el proceso ocupó más de 1 hora de tiempo de cómputo. Esta alternativa es recomendada en si resultase imprescindible conservar las propiedades topográficas en el DEM, pero debe tenerse en cuenta que su rendimiento es muy bajo. En los casos donde se use un DEM pequeño, se recomienda usar esta alternativa. Sin embargo, en esta investigación, con un DEM de más 300 millones de celdas, preferimos evaluar otras alternativas.
# Limpiar red manualmente en QGIS
# Aplicar v.clean, en QGIS
# Tallar la red densa
# time r.carve --overwrite --verbose raster=dem_pseudo_ortometrico \
# vector=red_mtn50k_cleaned output=dem_tallado depth=100
# echo "r.carve finalizado" | mail -s "r.carve finalizado" zoneminderjr@gmail.com
## real 320m40.373s # NO ALCANZÓ A TERMINAR USANDO RED DOMINICANA DEL MTN50K DENSA (red_mtn50k_cleaned)
# Tallar red de cursos largos
# time r.carve --overwrite --verbose raster=dem_pseudo_ortometrico \
# vector=red_mtn50k_cleaned_largos output=dem_tallado depth=100
# echo "r.carve finalizado" | mail -s "r.carve finalizado" zoneminderjr@gmail.com
## real 97m3.970s # COMPLETADO: USANDO RED DE CURSOS LARGOS SOLAMENTE (red_mtn50k_cleaned_largos)Alternativa 2. Tallado con GRASS GIS, usando álgebra de mapas r.mapcalc (elegida)
Tallado eficiente usando álgebra de mapas. Normalizamos el DEM, generamos una capa booleana ráster con la red de cursos largos, la restamos al DEM normalizado y luego multiplicamos el ráster resultante de la resta nuevamente por el rango del DEM (máximo - mínimo) para restablecer los valores originales fuera de las áreas talladas. El resultado es un DEM tallado, en el que los píxeles tallados (por donde circula la red) tenían una profundidad (valor negativo) equivalente al rango.
# Implementar esta alternativa con mapcalc de GRASS GIS y evaluar rendimiento
# https://www.youtube.com/watch?v=jHT_StPb_oM
# Limpiar red manualmente en QGIS
# Aplicar v.clean, en QGIS
# Tallar
# Rasterizar red (los píxeles de la red valdrán 1, el resto, nulo)
v.to.rast --overwrite input=red_mtn50k_cleaned_largos type=line use=val output=red_mtn50k_cleaned_largos
# Convertir nulos a cero
r.null map=red_mtn50k_cleaned_largos null=0
# Determinar estadísticas univariantes del DEM
r.univar map=dem_pseudo_ortometrico
# minimum: -51.4456
# maximum: 3102.34
# Aplicar normalización y resta
r.mapcalc --overwrite << EOF
eval(stddem = (dem_pseudo_ortometrico - -51.4456) / (3102.34 - -51.4456), \
stddemburn = stddem - red_mtn50k_cleaned_largos)
dem_tallado = (stddemburn * (3102.34 - -51.4456)) - 51.4456
# dem_tallado = stddemburn * dem_pseudo_ortometrico # Alternativa
EOF
echo "Tallado finalizado" | mail -s "Mensaje sobre tallado" zoneminderjr@gmail.com
FIGURA 8: DEM ALOS PALSAR sin aplicación de hidrografía (A), y con aplicación de hidrografía seleccionada (B). El DEM se representa como relieve sombreado y la aplicación se denota como un grabado oscurecido (cañón del río Payabo, Los Haitises, y río Yuna (proximidades de Arenoso, nordeste de República Dominicana)
Alternativa 3. Tallado con WhiteboxTools (descartada)
Como alternativa de procesamiento 3 usamos la herramienta FillBurn, basada en Saunders (2000) e implementada por John B. Lindsay (2016) en de WBT. Esta realiza dos modificaciones a la vez sobre el DEM; por una parte, graba la red, usando una profundidad por defecto y, por otro, rellena las depresiones. La herramienta mostró mejor rendimiento que la de GRASS GIS en cuanto a tiempo de cómputo, tanto con la red densa, que sí pudo grabarla sobre el DEM, como con la de cursos largos.
En particular, sobre el grabado de la red densa, notamos que los trazados no se alineaban bien con el DEM por problemas de alineación y ajuste de datums entre el mapa topográfico nacional y el DEM (NAD27 y WGS84). Como consecuencia, la red se grababa en topografías claramente muy accidentadas o a contrapendiente, generando así una red con duplicados y falsos positivos. No obstante, intentamos generar una red con dicho DEM, pero el resultado no fue satisfactorio en términos hidrológicos, por lo que optamos por probar con la red de cursos largos para generar el DEM tallado. En este caso, el DEM resultante fue muy diferente al original, especialmente en las áreas con depresiones. Por esta razón, decidimos usar el DEM tallado con la red de cursos largos generado por GRASS GIS.
# Alternativa 2: rápida, pero produce un tallado muy profundo y rellena depresiones
# Exportar dem_pseudo_ortometrico a GTiff con compresión LZW
# time r.out.gdal --overwrite --verbose createopt="COMPRESS=LZW,BIGTIFF=YES" \
# input=dem_pseudo_ortometrico \
# format=GTiff type=Float64 output=dem_pseudo_ortometrico.tif
# echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 1m0.248s
# Exportar red_mtn50k_cleaned.gpkg a shapefile
# ogr2ogr(src_datasource_name = '/media/jose/datos/alos-palsar-dem-rd/dem/red_mtn50k_cleaned.gpkg',
# dst_datasource_name = '/media/jose/datos/alos-palsar-dem-rd/dem/red_mtn50k_cleaned.shp',
# verbose=TRUE)
# Tallar con WBT
# time ~/WhiteboxTools_linux_amd64/WBT/whitebox_tools \
# --wd='/media/jose/datos/alos-palsar-dem-rd/dem/' \
# --run=FillBurn --dem='dem_pseudo_ortometrico.tif' \
# --streams=red_mtn50k_cleaned.shp --output='dem_tallado.tif' -v
# echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 9m21.980s
# Importar a GRASS GIS
# time r.import --overwrite input=dem_tallado.tif output=dem_tallado
# echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 0m38.519sDepresiones
Para producir límites de cuencas y redes de drenaje con sentido, los algoritmos de análisis hidrológico requieren que se identifiquen las depresiones capaces de capturar la escorrentía superficial (e.g. ponor, pérdidas). Usando la capa de litologías de República Dominicana (Mollat et al. 2004), identificamos y separamos las calizas que presentaban suficiente grado de karstificación, de acuerdo con nuestra experiencia de campo. Posteriormente, generamos una capa de depresiones con el complemento r.geomorphon usando el DEM como insumo (Jarosław Jasiewicz y Stepinski 2013). Adicionalmente, digitalizamos algunas depresiones cuya localización conocíamos por experiencia de terreno. Finalmente, intersectamos las tres fuentes para producir una capa comprensiva de las depresiones que capturan el flujo superficial.
FIGURA 9: “Geomórfonos” de República Dominicana generados a partir de DEM ALOS PALSAR. En cartela, detalle del cañón del río Payabo
No obstante, nuestro resultado debe tomarse con cautela en el relieve kárstico. Como bien es sabido, no todas las calizas representadas en la geología dominicana están lo suficientemente karstificadas como para desarrollar depresiones. Por esta razón, usamos la capa de calizas a discreción, y sólo conservamos aquellos afloramientos de calizas donde no se evidenciaba escorrentía superficial, y donde encontramos evidencia de depresiones en la topografía detallada y en imágenes satelitales. No obstante, gran parte de este trabajo se realizó manualmente, por lo que la exhaustividad no está del todo garantizada. Además, es virtualmente imposible identificar todas las depresiones que funcionan como pérdidas en imágenes satelitales o en mapas topográficos. Finalmente, un elemento adicional complica aún más las cosas en los relieves kársticos: muchas pérdidas no ocurren a través de una depresión topográficamente visible, pues gran parte de la infiltración se produce a través de fracturas en la roca, pasando al endokarst y a la zona vadosa de manera “silenciosa”, sin que veamos desde el aire la típica morfología deprimida (e.g. dolina).
# Crear geomórfonos
# WBT
# time ~/WhiteboxTools_linux_amd64/WBT/whitebox_tools \
# -r=Geomorphons -v --wd='/media/jose/datos/alos-palsar-dem-rd/dem/' \
# --dem=dem_tallado.tif -o=geomorfonos.tif --search=25 \
# --threshold=0 --tdist=0.0 --forms
# echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 6m52.298s #MUY EFICIENTE
# GRASS GIS
# time r.geomorphon --overwrite --verbose elevation=dem_tallado forms=geomorfonos search=25
time r.geomorphon --overwrite --verbose elevation=dem_pseudo_ortometrico forms=geomorfonos search=25
echo "r.geomorphon finalizado" | mail -s "Mensaje sobre r.geomorphon" zoneminderjr@gmail.com
## real 33m16.508s #MUY LENTO
# Extraer depresiones desde geomorfonos
r.mapcalc --overwrite expression="'depresiones_geomorfonos' = if(geomorfonos == 10, 1, null())"
# Importar depresiones manualmente digitalizadas a base de datos de GRASS GIS
v.import --overwrite input=depresiones_digitalizadas.gpkg output=depresiones_digitalizadas
# Convertir depresiones digitalizadas manualmente a ráster
v.to.rast --overwrite input=depresiones_digitalizadas type=area use=val output=depresiones_digitalizadas
# Importa la capa de calizas con depresiones en RD (de Mapa Geológico 250K)
v.import --overwrite input=calizas_con_depresiones.gpkg output=calizas_con_depresiones
# Convertir la capa de calizas con depresiones a ráster
v.to.rast --overwrite input=calizas_con_depresiones type=area \
use=val output=calizas_con_depresiones
# Adjuntar depresiones digitalizadas manualmente con caliz
r.mapcalc --overwrite \
expression="'depresiones_geomorfonos_calizas' = 'depresiones_geomorfonos' * 'calizas_con_depresiones'"
# Unir todas las depresiones en un único mapa
r.patch --overwrite input=depresiones_geomorfonos_calizas,depresiones_digitalizadas \
output=depresiones_todas
FIGURA 10: DEM ALOS PALSAR representado como mapa hipsómétrico (rojo y marrón representan terreno elevado, verde y azul claro terreno bajo) sobre relieve sombreado, mostrando el área de Guaraguao, Los Haitises, al sur del río Yuna (nordeste de República Dominicana). (A) sin mostrar depresiones, (B) mostrando depresiones en tonalidad azul oscuro
Extracción de red de drenaje y cuencas hidrográficas: análisis hidrológicos y morfométricos, jerarquía de red
El dilema: r.watershed o r.stream*
El complemento r.watershed es uno de los más completos de GRASS GIS para realizar análisis hidrológico a nivel de cuenca e incluso para redes. Con este complemento se pueden obtener los mapas de acumulación de flujo, talwegs, redes de drenaje y cuencas, entre otros. Estos productos son usados tanto en el análisis de cuenca, como en el análisis de red de drenaje (e.g. análisis hortoniano). Sin embargo, otro complemento de GRASS GIS se especializa en la extracción y caracterización de redes de drenaje. La elección de uno u otro dependerá de nuestro objetivo. Si nos interesa analizar cuencas hidrográficas, entonces, el complemento a usar será r.watershed. Si nos interesa analizar redes de drenaje y jerarquía hidrográfica, entonces tenemos dos opciones: r.watershed y r.stream*, ya sea combinados o de forma excluyente.
La familia r.stream* incluye un generador de red (r.stream.extract), un calculador de jerarquía de red (r.stream.order) y un generador de cuencas hidrográficas consecuente con la red de drenaje (r.stream.basins) (JarosŁaw Jasiewicz y Metz 2011). El complemento r.stream.order es bastante potente, porque nos permite caracterizar la red hidrográfica de manera muy completa usando la topología para determinar la jerarquía ordenada. Para correr r.stream.order, necesitamos el mapa de la red propiamente (stream_rast), y el de dirección de drenaje (direction). Por otra parte, r.stream.basins utliza como entrada los mapas producidos por r.stream.order para delimitar cuencas de acuerdo con la ordenación jerárquica elegida.
Tanto r.stream.extract como r.waterhsed son capaces de generar productos que sirven de insumo a otros complementos, pero con ligeras diferencias. En un uso combinado de estos complementos, para garantizar la consistencia del resultado, la clave está en no combinar un producto de r.watershed con uno de r.stream.extract. Por ejemplo, para calcular órdenes de red con r.stream.order, no se recomienda combinar el mapa de la red producido por r.watershed con el de dirección de drenaje de r.stream.extract. Ambas entradas deben ser producidas por un mismo algoritmo.
Por lo tanto, en este punto tenemos dos opciones:
Opción 1: generar tanto
stream_rastydirectionconr.watershed.Opción 2: generar tanto
stream_rastydirectionconr.stream.extract.
Dado que nuestro interés principal es la jerarquía de red, pues la utilizamos como criterio de selección de sitios idóneos para estaciones hidrométricas, podríamos saltar directamente a ejecutar el complemento r.stream.extract, y generar productos que usaremos posteriormente en r.stream.order. Esto nos ahorraría el paso de ejecutar r.watershed. Ahora bien, dado que nos interesa que r.stream.order calcule la jerarquía por el método Horton, entonces necesitaremos el mapa de acumulación, el cual sólo puede producirlo r.watershed. Comencemos precisamente con ello, generando el mapa de acumulación.
Dado que a partir de este bloque de código, inician los análisis hidrológicos, antes de ejecutar r.watershed, aplicamos una máscara ajustada a la línea de costa y los límites fronterizos del país. Esta máscara más ajustada asegura que las redes extraídas no se extiendan hacia en el mar. Por otra parte, para evitar interrupciones abruptas y no realistas del flujo, creamos una zona de influencia en el límite fronterizo para permitir la salida y entrada de flujo a través de éste.
# Importar máscara
v.import --overwrite input=mascara-1km-solo-en-frontera.gpkg output=mascara_1km_solo_en_frontera
# Fijar máscara (EJECUTAR SÓLO SI ES ESTRICTAMENTE NECESARIO, PUES TARDA MUCHO)
r.mask -r
r.mask vector=mascara_1km_solo_en_fronteraIniciamos los análisis hidrológicos con la nueva máscara.
time r.watershed --overwrite --verbose elevation=dem_tallado \
depression=depresiones_todas accumulation=rwshed_acum \
threshold=180 stream=rwshed_talwegs \
# drainage=rwshed_direccion_drenaje basin=rwshed_cuencas half_basin=rwshed_hemicuenca \
# tci=rwshed_tci spi=rwshed_spi \
# length_slope=rwshed_longitud_vertiente slope_steepness=rwshed_empinamiento \
# retention=rwshed_retencion_flujo max_slope_length=rwshed_max_longitud_vertiente \
memory=32000
echo "r.watershed finalizado" | mail -s "Mensaje sobre r.watershed" zoneminderjr@gmail.com
## real 9m47.041s
FIGURA 11: Mapas de acumulación de flujo y red extraída con r.watershed
Extracción de red de drenaje y cuencas hidrográficas según orden jerárquico
Con el DEM y la capa de acumulación generada por r.watershed, extrajimos la red hidrográfica, para lo cual, previamente elegimos, por medio de inspección visual, un umbral óptimo de acumulación. La elección del umbral óptimo puede variar en función de las características específicas de cada cuenca, como la pendiente y tamaño de la misma, la litología, el clima, entre otros atributos. En consecuencia, la práctica habitual recomienda realizar varias corridas de análisis con diferentes valores umbral para obtener la mejor configuración de la red hidrográfica para el área de interés. Los criterios que fijamos para elegir un umbral óptimo de acumulación fueron los siguientes:
Umbrales consistentes con estudios similares.
Densidad de red suficiente como para producir una red representativa.
Evitar la generalización de la red al punto que invisibilice lugares idóneos.
Evitar una red extremadamente densa con lugares que no reúnen las mínimas características hidrológicas (e.g. cursos con caudal permamennte o semipermanente).
Dado que nuestro DEM cuenta con una resolución espacial de 12.5 m, exploramos diferentes umbrales para determinar la extensión de las cuencas hidrográficas y la densidad de la red hidrográfica que nos resultara idónea. Seleccionamos los valores 180, 540 y 900 celdas como umbrales de acumulación en el complemento r.stream.extract, los cuales equivalen, respectivamente, a 3, 8 y 14 hectáreas en términos de superficie, respectivamente. Estos umbrales se inscriben en el rango de valores usados en estudios consultados (Marchesini et al. 2021), donde evalúan por métodos de “áreas susceptibles a inundación” (Freedman y Diaconis 1981), el rango de umbrales a elegir.
# Extraer redes de drenaje para tres umbrales de acumulación distintos
# En bucle
for i in `echo {180..900..360}`; \
do echo -e "\nTRABAJANDO EL UMBRAL DE ACUMULACIÓN $i ...\n"; \
time r.stream.extract --overwrite elevation=dem_tallado accumulation=rwshed_acum \
depression=depresiones_todas threshold=$i \
stream_vector=rstream_talwegs_umbral_$i stream_raster=rstream_talwegs_umbral_$i \
direction=rstream_direccion_umbral_$i memory=32000; \
echo -e "r.stream.extract umbral $i finalizado" | mail -s "Mensaje sobre r.stream.extract" zoneminderjr@gmail.com; \
done
## real 11m28.455s
## real 11m26.908s
## real 11m30.074s
# Único umbral, para testing
# time r.stream.extract --overwrite elevation=dem_tallado accumulation=rwshed_acum \
# depression=depresiones_todas threshold=64 \
# stream_vector=rstream_talwegs_umbral_64 stream_raster=rstream_talwegs_umbral_64 \
# direction=rstream_direccion_umbral_64 memory=32000
# echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 11m46.930s
# Calcular estadisticos, y pasar a archivo
for i in `echo {180..900..360}`; \
do v.to.db -c -p option=length map=rstream_talwegs_umbral_$i | grep total > stats_length_rstream_talwegs_umbral_$i.txt; \
done
for i in `echo {180..900..360}`; \
do v.to.db -c -p option=count map=rstream_talwegs_umbral_$i | grep total > stats_count_rstream_talwegs_umbral_$i.txt; \
doneA partir del archivo de texto generado, obtuvimos los estadísticos básicos: se trata de una red compuesta por 855503, 391841 y 265491 segmentos, respectivamente para cada umbral de 180, 540 y 900 celdas, los cuales suman totales de 138444.33, 98213.36 y 82080.67 kilómetros de longitud, respectivamente.
Tras realizar una inspección visual de los resultados producidos por el complemento r.stream.extract, elegimos los resultados correspondientes al umbral de 540 celdas (8 ha), por adaptarse mejor a los fines de nuestro estudio. Las redes producidas con este umbral resultaron ser bastante limpias y diversas en el territorio dominicano, mostrando una hidrografía representativa sin ocultar patrones de variabilidad, necesarios en la toma de decisiones sobre selección de sitios para instalar estaciones hidrométricas. Además, como se muestra a continuación, la red generada con este umbral, alcanzó un orden máximo 8 según los métodos de Strahler y Horton, lo cual nos permitió usar, de manera preferente,pero no exclusiva, los órdenes de red 5 y 6 como candidatos idóneos para el establecimeinto de estaciones hidrométricas, por su mayor probabilidad de ser ríos permanentes o semipermanentes, sin llegar a tener valles de fondo muy ancho.
FIGURA 12: Red de drenaje extraída para tres umbrales de acumulación: (A) 180 celdas, equivalente a ~3 ha; (B) 540 celdas, equivalente a ~8 ha; (C) 900 celdas, equivalente a ~14 ha. La imagen de fondo es un relieve sombreado a partir de DEM ALOS PALSAR, mostrando el área de El Arroyazo en la reserva científica Ébano Verde (provincia La Vega, cordillera Central de República Dominicana)
Generar órdenes de red
# Extraer orden de red
# En bucle
for i in `echo {180..900..360}`; \
do echo -e "\nTRABAJANDO EL UMBRAL DE ACUMULACIÓN $i ...\n"; \
time r.stream.order --overwrite stream_rast=rstream_talwegs_umbral_$i \
direction=rstream_direccion_umbral_$i \
elevation=dem_tallado accumulation=rwshed_acum \
stream_vect=rstream_orden_de_red_umbral_$i \
strahler=rstream_orden_strahler_de_red_umbral_$i \
horton=rstream_orden_horton_de_red_umbral_$i \
topo=topologia_orden_umbral_$i memory=32000; \
echo -e "r.stream.order umbral de acumulación $i finalizado" | mail -s "Mensaje sobre r.stream.order" zoneminderjr@gmail.com; \
done
## real 1m34.983s
## real 1m18.662s
## real 1m14.986s
# Único
# time r.stream.order --overwrite stream_rast=rstream_talwegs direction=rstream_direccion \
# elevation=dem_tallado accumulation=rwshed_acum stream_vect=order_todos \
# topo=topologia_orden memory=32000
# echo "Job finished" | mail -s "Job finished" zoneminderjr@gmail.com
## real 
FIGURA 13: Representación de la red hidrográfica dominicana usando simbolizando el grosor de los segmentos en función de su orden de red (método de Strahler). El orden máximo alcanzado fue de 8. Esta red fue generada usando un umbral de acumulación de 540 celdas (~8 ha)
FIGURA 14: Orden de red de Strahler para redes de drenaje generadas a partir de tres umbrales de acumulación: (A) 180 celdas, equivalente a ~3 ha; (B) 540 celdas, equivalente a ~8 ha; (C) 900 celdas, equivalente a ~14 ha. El área mostrada corresponde al río San Juan, afluente del río Yaque del Sur (vertiente sur de la cordillera Central de República Dominicana)
FIGURA 15: Orden de red de Strahler en el área del pico de la Viuda y Sabana Vieja, provincia San Juan (vertiente sur de la cordillera Central de República Dominicana). Esta red fue generada usando un umbral de acumulación de 180 celdas, equivalente a ~3 ha. De fondo, mapa topográfico nacional escala 1:50,000 y relieve sombreado
Delimitar cuencas y subcuencas según la jerarquía de red
Usamos el complemento r.stream.basins dentro de un bucle for para delimitar las cuencas y subcuencas según su orden de red máximo. Las unidades delimitadas por este procedimiento incluyen tanto cuencas y subcuencas, por lo que, la mayoría contiene redes de drenaje tributarias de otro río.
# Delimitar cuencas según jerarquía
# En bucle
for i in `echo {180..900..360}`; \
do for j in `echo {1..8..1}`; \
do echo -e "\nTRABAJANDO EL UMBRAL DE ACUMULACIÓN $i, orden $j...\n"; \
time r.stream.basins -c --overwrite direction=rstream_direccion_umbral_$i \
stream_rast=rstream_orden_strahler_de_red_umbral_$i cats=$j \
basins=rstream_cuencas_strahler_umbral_${i}_orden_$j memory=32000; \
done; \
echo -e "r.stream.basins umbral de acumulación $i finalizado" | mail -s "Mensaje sobre r.stream.basins" zoneminderjr@gmail.com; \
done
## real ~ 0m40s repetido tantas veces como órdenes para cada umbral de acumulaciónDelimitar cuencas con drenaje final
En esta sección aplicamos el mismo complemento, pero esta vez para delimitar las cuencas con drenaje final, es decir, cuencas exorreicas (drenan al mar) o endorreicas (desembocan en lagos, lagunas o en pérdidas del karst, como ponors). Es decir, se trata de cuencas propiamente en la acepción más formal del término, que significa que no existe (no se conoce) prolongación del drenaje fuera de ellas.
# Delimitar cuencas terminales
# En bucle
for i in `echo {180..900..360}`; \
do for j in `echo {1..8..1}`; \
do echo -e "\nTRABAJANDO EL UMBRAL DE ACUMULACIÓN $i, orden $j...\n"; \
time r.stream.basins -lc --overwrite direction=rstream_direccion_umbral_$i \
stream_rast=rstream_orden_strahler_de_red_umbral_$i cats=$j \
basins=rstream_cuencas_strahler_terminal_umbral_${i}_orden_$j memory=32000; \
done; \
echo -e "r.stream.basins umbral de acumulación $i finalizado" | mail -s "Mensaje sobre r.stream.basins" zoneminderjr@gmail.com; \
done
## real ~ 0m40s repetido tantas veces como órdenes para cada umbral de acumulaciónConvertir cuencas a vectorial. Sólo las de umbral 540
En esta sección, convertimos las cuencas anteriores a vectorial para su mejor manejo y representación. Asimismo, representamos el mapa de las estaciones hidrométricas existentes sobre las cuencas categorizadas en función del orden de red máximo encontrado en ellas. A primera vista, notamos que la cobertura de la red es relativamente baja y que, a priori, existe un déficit importante de estaciones hidrométricas en la mayoría de las cuencas de orden 5 y superiores, con una marcada concentración en cuencas con presas. Reconocemos que el establecimiento de estaciones hidrométricas asociadas a las presas, es una labor de rutina y necesaria, pero el desbalance entre cuencas instrumentadas y no instrumentadas es significativo, lo cual en última instancia contribuye al sesgo en el dato hidrométrico. El mapa muestra tanto las estaciones que se encontraban en estado “Regular” como la de estado “Bueno” del informe de INDRHI de 2019 (INDRHI 2019). Asimismo, añadimos al mapa las dos estaciones manejadas por el Servicio Geológico Nacional (SGN).
for i in `echo {1..8..1}`; \
do r.to.vect --overwrite input=rstream_cuencas_strahler_terminal_umbral_540_orden_$i \
output=rstream_cuencas_strahler_terminal_umbral_540_orden_$i type=area; \
v.db.addcolumn rstream_cuencas_strahler_terminal_umbral_540_orden_$i columns="strahler int"; \
v.db.update rstream_cuencas_strahler_terminal_umbral_540_orden_$i \
col=strahler value=$i where="strahler IS NULL"; \
done
v.patch -e --overwrite \
input=`g.list type=v pattern='rstream_cuencas_strahler_terminal_umbral_540_orden_*' separator=comma` \
output=rstream_cuencas_strahler_terminal_umbral_540_todos
FIGURA 16: Cuencas de orden 4 o superior de República Dominicana, extraídas a partir de una red de drenaje creada usando 540 celdas (aprox. 8 ha) como umbral de acumulación del DEM ALOS PALSAR. Nótese el marcado déficit de estaciones en gran parte del territorio
Áreas de prospección de sitios idóneos para estaciones hidrométricas
Las redes de drenaje y las cuencas hidrográficas generadas en los pasos anteriores, constituyen la base primordial del presente estudio y uno de sus principales objetivos. Con la red generada, definimos el área de interés que nos permitió seleccionar sitios para estaciones hidrométricas, lo cual constituye el segundo objetivo del presente estudio. Explicamos en esta sección cómo desarrollamos ese segundo objetivo.
Con la red disponible, generando algunas variables adicionales, podríamos aplicar, de forma un tanto ingenua, algoritmos de reclasificación de variables y álgebra de mapas para producir una lista sitios que cumplan los criterios fijados.
Con la red disponible, y generando algunas variables adicionales, podríamos aplicar algoritmos de reclasificación de variables y álgebra de mapas para producir una lista de sitios que cumplan con los criterios establecidos. Sin embargo, esta aproximación podría considerarse un tanto ingenua y hasta ineficiente, pues estaríamos ignorando los resultados que obtuvimos en el estudio complementario a éste titulado “Selección de sitios para el establecimiento de una red de estaciones meteoclimáticas en República Dominicana usando decisión multicriterio y análisis de vecindad”.
Para focalizar la búsqueda, nuestra estrategia consistió en generar lo que denominamos “áreas de prospección”. Estas áreas son porciones manejables del territorio que rodean a estaciones meteoclimáticas existentes y propuestas por nosotros en el mencionado estudio (e.g. zonas buffer). Así, dentro de las áreas de prospección, generamos múltiples variables y aplicamos criterios de manera eficiente, al focalizadr en esos pequeños espacios del territorio dominicano nuestros recursos de hardware y de software.
Esta estrategia tuvo una doble ventaja. Por un lado, redujo los tiempos de cómputo y, por otro, al restringir la búsqueda al entorno de un sitio con una estación meteoclimática (existente o propuesta), aprovechamos la sinergia de la co-observación del dato climático e hidrométrico. Adicionalmente, esto nos permitió circunscribir la búsqueda a áreas que ya fueron priorizadas por otros criterios, además de los meramente hidrográficos. En el referido estudio sobre estaciones meteoclimáticas, consolidamos una relación exhaustiva de puntos con coordenadas para distintos escenarios de densidad de estaciones. En el caso de las estaciones hidrométricas, dichos escenarios también resultan convenientes, pues ofrecen a las tomadores de decisiones distintas alternativas de inversión.
Materializamos nuestras áreas prospección como círculos de radio \(R_p\), el cual debemos obtener, que cubren los sitios de las estaciones meteoclimáticas ya existentes y propuestas. El estudio garantizaba que los sitios se espaciaran adecuadamente entre sí, para garantizar cubrir un área \(A\) representativa (normalmente, una superficie sugerida por la OMM), y evitando al mismo tiempo redundancia entre estaciones. Si asumimos que la geometría de áreas de representación es circular, entonces la separación promedio entre sitios de estaciones climáticas equivale al diámetro \(D\) de un círculo imaginario con área \(A\).
Dado que nuestro objetivo es localizar sitios idóneos para establecer estaciones hidrométricas, decidimos que nuestras áreas de prospección fuesen círculos de \(R_p=D/3\) centrados en las estaciones meteoclimáticas propuestas y existentes. Con esto flexibilizamos la búsqueda, permitiéndonos seleccionar sitios dentro de un área lo suficientemente grande, sin comprometer a su vez el criterio de separación entre estaciones para evitar redundancia.
Un ejemplo práctico ilustra el cálculo de \(R_p\). Usemos el escenario “una estación por cada 250 km2”. Esto significa que \(A=\) 250 km2, entonces el diámetro \(D\) se obtiene por: \(D=\frac{\sqrt{4A/pi}}{3}=\frac{\sqrt{4 \times 250/3.1416}}{3}\approx 6\) km. Para otros escenarios, e.g. una estación por cada 150 km2, bastaría con sustituir \(A\) por 150 en la fórmula.
Generamos las áreas de prospección específicamente con el siguiente código, las cuales pueden verse representadas en el mapa siguiente. Con el bloque importamos primero la propuesta de sitios de estaciones meteoclimáticas para cada escenario de densidad (100, 150 y 250 km2), considerando sólo la red donde se eliminó redundancia por proximidad con estaciones meteoclimáticas de INDRHI y ONAMET existentes que se encontraban en estado “bueno”. Posteriormente, importamos la red de estaciones meteoclimáticas de INDRHI y ONAMET existentes, conservando sólo aquellas en estado bueno. Finalmente, combinamos, en un único vectorial, las estaciones y los sitios propuestos, para generar un buffer que serán nuestas áreas de prospección.
# Propuestas de estaciones
vareaskm2=(100 150 250)
for areakm2 in "${vareaskm2[@]}"; do v.import --overwrite \
input=escenario_${areakm2}_km2_por_estacion_exclusion_redundancia_activas_buenas.gpkg \
output=escenario_${areakm2}_km2;
done
# Estaciones climáticas del INDRHI
v.import --overwrite input=con_indicacion_estatus_climaticas_indrhi.gpkg \
output=tmp; \
v.extract --overwrite input=tmp where='Estado=="Bueno"' output=estaciones_indrhi
g.remove -f type=vector name=tmp
# Estaciones climáticas de ONAMET
v.import --overwrite input=con_indicacion_estatus_onamet.gpkg \
output=tmp; \
v.extract --overwrite input=tmp where='estado=="activa"' output=estaciones_onamet
g.remove -f type=vector name=tmp
#Combinar estaciones y sitios propuestos en un único vectorial
for areakm2 in "${vareaskm2[@]}"; do v.patch --overwrite \
input=escenario_${areakm2}_km2,estaciones_indrhi,estaciones_onamet \
output=puntos_para_areas_prospeccion_${areakm2}_km2;
done
# Crear buffer de tamaño variable, con categorías únicas ("cat" secuencial)
for areakm2 in "${vareaskm2[@]}"; do v.buffer --overwrite \
input=puntos_para_areas_prospeccion_${areakm2}_km2 distance=`echo "scale=2;(sqrt(4*$areakm2/3.14159))*1000/3" | bc` \
output=areas_prospeccion_${areakm2}_km2;
v.category --overwrite input=areas_prospeccion_${areakm2}_km2 output=tmp option=del cat=-1;
v.category --overwrite input=tmp output=areas_prospeccion_${areakm2}_km2 option=add step=1;
v.db.addtable areas_prospeccion_${areakm2}_km2 columns="cat int";
doneLos pasos siguientes consistieron en generar, dentro de las áreas de prospección, los diferentes criterios que aplicamos a la selección de sitios para el establecimiento de estaciones hidrométricas Para ello, ejecutamos los complementos correspondientes de GRASS GIS (e.g. r.valley.bottom) fijando como máscara las áreas de prospección creadas arriba. Este paso fue crucial, pues conseguimos eficientizar la generación de las variables de entrada para seleccionar sitios prioritarios.
Generación de criterios
Elegimos los siguientes criterios de priorización de sitios para el establecimiento de estaciones hidrométricas, los cuales calculamos exclusivamente dentro de las áreas prospección:
- Orden de red: 4 para órdenes 4, 5, 6; 1 para órdenes ≤ 3 o ≥ 7.
- Criterios de Rantz seleccionados:
- Fondo de valle: 4 para fondos de valle pequeños; 1 para fondos de valle grandes.
- Sinuosidad del tramo (#1 de lista de Rantz): 4 para rectilíneo (sinuosidad baja); 1 para curvo (sinuosidad alta)
- Tipo de sustrato rocoso, rocas consistentes (#5): 4 para magmáticas nada o poco alteradas; 3 para metamórficas (no mármoles); 2 para sedimentarias (no calizas); 1 para calizas.
- Confluencias y efectos mareales (#7): 4 para distante; 1 para próximo
- Accesible (#9): 4 para próximo; 1 para distante.
- Instrumentación en la subcuenca específica: se aplica al finalizar, realizando un análisis de vecindad por subcuenca, para determinar si el sitio propuesto se encuentra próximo a otro existente.
Fondo de valle
Usando el complemento `r.valley.bottom, obtuvimos el “índice múltiresolución de planitud de fondo de valle” (MrVBF), con el cual estimamos la superficie ocupada por este. Este criterio nos sirvió para discriminar entre talwegs con fondo de valle anchos, (menos prioritarios, por su mayor inversión requerida), y talwegs con fondos de valle estrechos (prioritarios).
Nota sobre la eficiencia del complemento
r.valley.bottom. En un primer intento, ejecutamos el siguiente bloque de código para obtener el MrVBF de todo el país. Entendíamos posible lograrlo y, una vez obtenido, recortarlo en las áreas de prospección, un procedimiento que hemos realizado con muchos otros algoritmos, pero al parecer, el complementor.valley.bottomno está preparado para demandas tan exigentes. Nos vimos en la necesidad de buscar formas distintas de cómputo, porque la ejecución del algoritmo para todo el país tardó más de 4 horas, y el ráster resultante sólo contenía valores sin datos, lo cual atribuimos a la necesaria definición de parámetros específicos para un territorio topográficamente muy diverso. Esto nos obligó a calcular un ráster de MrVBF para sectores más pequeños, en concreto, sólo para nuestras áreas de prospección.
# Tarda un tiempo considerable
# Al elegir min_cells=10000 (por omisión asumiría 2), el número de pasos es ocho,
# que es mucho menor que los 12 o 16 que ejecutaría con valores por omisión.
# Sin embargo, la última prueba exitosa, consistió en utilizar, como valor
# del parámetro min_cells, el número de celdas de la región, en este caso 911808546
# Esto redujo el número de pasos a 3, y el resultado fue satisfactorio.
vareaskm2=(100 150 250)
for areakm2 in "${vareaskm2[@]}"; do \
r.mask -r; \
r.mask vector=areas_prospeccion_${areakm2}_km2; \
time r.valley.bottom --verbose --overwrite elevation=dem_pseudo_ortometrico \
mrvbf=rvb_fondo_de_valle_${areakm2}_km2 \
mrrtf=rvb_crestas_${areakm2}_km2 \
min_cells=911808546; \
echo -e "Fondo de valle finalizado" | mail -s "Mensaje sobre r.valley.bottom"; \
echo -e "Fondo de valle para $i km2 finalizado" | mail -s "Mensaje sobre r.valley.bottom" zoneminderjr@gmail.com; \
done
# real 82m3.368s
# real 81m34.817s
# real 80m59.169sNota sobre criterios no usados en esta primera evaluación, reservados para futuros intentos
Transcribimos en esta sección una lista de criterios inicialmente considerados, pero que finalmente decidimos no incluir en esta primera evaluación, a efectos de mantener el procedimiento lo más simple posible, y debido a que el tiempo de cómputo que requerían excedía el inicialmente previsto.
- Frecuencia de inundado. Se puede estimar por teledetección óptica y radar, pero consideramos que no aportaría nuevos elementos más de allá de la propia inspección visual que realizamos en imágenes de alta resolución. Para más detalles, consultar el trabajo de Chen y Zhao (2022), y ver repositorio asociado.
- SPI, ¿balance hídrico? Esta variable se podría generar para el país, pero es probable que no aporte muchos más elementos para elegir sitios prioritarios. De todas formas, se podría evaluar su cómputo.
- Variables que pueden generarse con el complemento
r.basiny otros de GRASS GIS y software que hemos desarrollados en otras investigaciones (Martinez Batlle 2019; Martínez-Batlle 2019). Esta variables son relativamente fáciles de generar, pero abultarían mucho la selección de criterios. No obstante, para las cuencas drenadas por las estaciones hidrométricas actuales y propuestas, podría ser de gran interés como forma de caracterizar el territorio instrumentado y potencialmente a instrumentar:- Tiempo de concentración (alejamiento medio).
- Coeficiente de torrencialidad.
- Razón de elongación / compacidad / circularidad.
- Concavidad transversal.
- TPI, TRI, TWI.
- Pendiente media.
- Índice de concavidad.
- Integral hipsométrica.